stanEx5
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(4,1),mar=c(1,5,1,1))
drw[,1,] |> plot(type='l',xlab='',ylab=pr)
drw[,2,] |> plot(type='l',xlab='',ylab=pr)
drw[,3,] |> plot(type='l',xlab='',ylab=pr)
drw[,4,] |> plot(type='l',xlab='',ylab=pr)
par(mar=c(3,5,3,3))
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}hierarchical model
class c=1~k
b0c~N(b00,sb0)
b1c~N(b10,sb1)
yc~N(b0c+b1c*x,s)n=50
k=9
x=runif(n,0,20)
c=sample(letters[1:k],n,replace=T)
b00=rnorm(k,10,5)
b0=b00[factor(c)]
b10=rnorm(k,2,1)
b1=b10[factor(c)]
y=rnorm(n,b0+b1*x,5)
qplot(x,y,shape=c,size=I(2))+
scale_shape_manual(values=1:k)
qplot(x,y,col=c)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcitionestimate as no class
ex8-0.stan
y~N(b00+b10*x,s)
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}mdl=cmdstan_model('./ex8-0.stan')
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
fn(mdl,data)## Initial log joint probability = -199299
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpFvED78/model-9918b9cc8.stan', line 14, column 2 to column 22)
## Error evaluating model log probability: Non-finite gradient.
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpFvED78/model-9918b9cc8.stan', line 14, column 2 to column 22)
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpFvED78/model-9918b9cc8.stan', line 14, column 2 to column 22)
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 38 -156.345 0.0262739 0.00171095 1 1 86
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -156.34
## b0 6.78
## b1 1.84
## s 13.83
## m0[1] 35.77
## m0[2] 29.30
## m0[3] 21.80
## m0[4] 8.04
## m0[5] 43.23
## m0[6] 37.14
## m0[7] 8.68
## m0[8] 32.52
## m0[9] 36.62
## m0[10] 36.72
## m0[11] 18.70
## m0[12] 42.85
## m0[13] 43.41
## m0[14] 25.13
## m0[15] 12.52
## m0[16] 18.77
## m0[17] 11.51
## m0[18] 31.43
## m0[19] 7.65
## m0[20] 32.54
## m0[21] 31.77
## m0[22] 28.54
## m0[23] 42.99
## m0[24] 29.98
## m0[25] 43.39
## m0[26] 33.41
## m0[27] 42.39
## m0[28] 41.15
## m0[29] 13.95
## m0[30] 8.30
## m0[31] 25.88
## m0[32] 39.08
## m0[33] 40.99
## m0[34] 14.54
## m0[35] 21.13
## m0[36] 9.55
## m0[37] 30.21
## m0[38] 23.15
## m0[39] 26.80
## m0[40] 7.31
## m0[41] 9.39
## m0[42] 38.97
## m0[43] 19.64
## m0[44] 22.56
## m0[45] 35.67
## m0[46] 34.55
## m0[47] 34.54
## m0[48] 34.14
## m0[49] 25.23
## m0[50] 35.19
## y0[1] 35.35
## y0[2] 23.41
## y0[3] 23.58
## y0[4] -2.21
## y0[5] 45.05
## y0[6] 29.95
## y0[7] -17.55
## y0[8] 17.57
## y0[9] 26.08
## y0[10] 33.38
## y0[11] 7.80
## y0[12] 36.52
## y0[13] 29.11
## y0[14] 5.03
## y0[15] 0.58
## y0[16] -8.45
## y0[17] 13.32
## y0[18] 23.72
## y0[19] 17.24
## y0[20] 41.35
## y0[21] 4.75
## y0[22] 21.90
## y0[23] 46.17
## y0[24] 48.29
## y0[25] 16.93
## y0[26] 13.42
## y0[27] 45.88
## y0[28] 31.25
## y0[29] 15.02
## y0[30] 15.62
## y0[31] 19.26
## y0[32] 35.38
## y0[33] 37.60
## y0[34] 24.29
## y0[35] 26.81
## y0[36] -15.65
## y0[37] 44.23
## y0[38] 25.64
## y0[39] 9.98
## y0[40] -13.88
## y0[41] -3.61
## y0[42] 45.82
## y0[43] 5.23
## y0[44] 19.56
## y0[45] 51.28
## y0[46] 55.60
## y0[47] 18.45
## y0[48] 38.34
## y0[49] 11.28
## y0[50] 44.00
## m1[1] 7.31
## m1[2] 11.32
## m1[3] 15.33
## m1[4] 19.34
## m1[5] 23.35
## m1[6] 27.36
## m1[7] 31.37
## m1[8] 35.39
## m1[9] 39.40
## m1[10] 43.41
## y1[1] 6.80
## y1[2] 22.67
## y1[3] 21.33
## y1[4] 47.35
## y1[5] 30.68
## y1[6] 15.24
## y1[7] 14.36
## y1[8] 21.25
## y1[9] 33.90
## y1[10] 65.01
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.2 seconds.
## Chain 4 finished in 0.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -155.23 -154.91 1.27 0.96 -157.63 -153.89 1.01 487 470
## b0 6.72 6.74 4.23 4.05 0.24 13.63 1.03 221 275
## b1 1.84 1.84 0.33 0.32 1.30 2.38 1.02 268 274
## s 14.45 14.34 1.51 1.49 12.14 17.12 1.00 2133 1398
## m0[1] 35.72 35.72 2.51 2.54 31.61 39.69 1.00 1542 1184
## m0[2] 29.25 29.21 2.05 2.03 25.94 32.52 1.00 1149 980
## m0[3] 21.75 21.76 2.26 2.24 18.17 25.26 1.01 314 364
## m0[4] 7.98 8.00 4.04 3.88 1.74 14.59 1.03 222 274
## m0[5] 43.18 43.16 3.48 3.43 37.45 48.86 1.00 717 955
## m0[6] 37.09 37.12 2.67 2.65 32.71 41.30 1.00 1439 1166
## m0[7] 8.62 8.66 3.94 3.77 2.48 15.06 1.03 222 277
## m0[8] 32.47 32.46 2.22 2.22 28.87 36.06 1.00 1599 1337
## m0[9] 36.57 36.60 2.61 2.61 32.30 40.70 1.00 1495 1181
## m0[10] 36.67 36.70 2.62 2.62 32.38 40.82 1.00 1486 1181
## m0[11] 18.65 18.66 2.56 2.54 14.58 22.73 1.02 261 315
## m0[12] 42.80 42.77 3.43 3.38 37.15 48.36 1.00 739 996
## m0[13] 43.36 43.32 3.51 3.45 37.58 49.07 1.00 708 936
## m0[14] 25.08 25.13 2.06 2.12 21.79 28.36 1.00 457 448
## m0[15] 12.47 12.51 3.36 3.23 7.17 17.83 1.03 227 290
## m0[16] 18.71 18.73 2.55 2.53 14.66 22.78 1.02 261 315
## m0[17] 11.46 11.52 3.51 3.35 5.97 17.09 1.03 225 290
## m0[18] 31.38 31.37 2.14 2.16 27.91 34.85 1.00 1560 1364
## m0[19] 7.60 7.62 4.10 3.92 1.30 14.31 1.03 221 275
## m0[20] 32.49 32.47 2.22 2.21 28.88 36.08 1.00 1600 1337
## m0[21] 31.71 31.70 2.16 2.18 28.24 35.20 1.00 1575 1367
## m0[22] 28.49 28.44 2.03 2.03 25.20 31.73 1.00 896 951
## m0[23] 42.94 42.91 3.45 3.40 37.26 48.56 1.00 730 955
## m0[24] 29.92 29.92 2.07 2.06 26.56 33.25 1.00 1365 1275
## m0[25] 43.34 43.31 3.51 3.45 37.56 49.05 1.00 709 936
## m0[26] 33.35 33.35 2.29 2.28 29.62 37.07 1.00 1606 1280
## m0[27] 42.34 42.30 3.36 3.32 36.79 47.79 1.00 769 965
## m0[28] 41.10 41.06 3.19 3.14 35.81 46.25 1.00 865 1018
## m0[29] 13.90 13.93 3.16 3.06 8.90 18.92 1.03 231 292
## m0[30] 8.24 8.29 4.00 3.86 2.05 14.78 1.03 222 274
## m0[31] 25.83 25.89 2.04 2.11 22.58 29.08 1.00 514 518
## m0[32] 39.03 39.01 2.91 2.91 34.17 43.62 1.00 1135 1040
## m0[33] 40.95 40.89 3.16 3.12 35.68 46.05 1.00 881 1018
## m0[34] 14.48 14.50 3.08 2.98 9.60 19.39 1.03 233 299
## m0[35] 21.07 21.08 2.32 2.29 17.36 24.71 1.01 299 342
## m0[36] 9.49 9.52 3.80 3.62 3.54 15.73 1.03 223 283
## m0[37] 30.16 30.16 2.08 2.08 26.80 33.53 1.00 1415 1258
## m0[38] 23.10 23.15 2.16 2.17 19.68 26.48 1.00 356 388
## m0[39] 26.74 26.77 2.02 2.07 23.51 29.99 1.00 605 651
## m0[40] 7.25 7.27 4.15 3.96 0.89 14.02 1.03 221 275
## m0[41] 9.34 9.35 3.83 3.63 3.34 15.61 1.03 223 283
## m0[42] 38.92 38.92 2.89 2.89 34.09 43.49 1.00 1153 1040
## m0[43] 19.58 19.60 2.46 2.44 15.65 23.46 1.02 273 334
## m0[44] 22.51 22.54 2.20 2.16 19.01 25.93 1.01 336 378
## m0[45] 35.62 35.62 2.50 2.52 31.53 39.58 1.00 1546 1184
## m0[46] 34.50 34.51 2.39 2.37 30.57 38.32 1.00 1585 1296
## m0[47] 34.49 34.50 2.39 2.37 30.56 38.31 1.00 1585 1296
## m0[48] 34.09 34.09 2.35 2.33 30.24 37.88 1.00 1595 1278
## m0[49] 25.18 25.23 2.06 2.11 21.90 28.46 1.00 464 449
## m0[50] 35.14 35.15 2.45 2.45 31.14 39.01 1.00 1563 1200
## y0[1] 34.96 35.27 14.93 14.19 9.77 59.35 1.00 2063 1972
## y0[2] 29.36 29.27 14.74 14.45 5.57 53.34 1.00 2055 1769
## y0[3] 21.39 21.26 14.31 14.01 -1.60 45.12 1.00 1840 1905
## y0[4] 7.72 7.74 14.69 14.78 -16.94 31.71 1.00 1789 1731
## y0[5] 43.17 43.42 14.59 14.04 18.97 66.35 1.00 1694 1704
## y0[6] 37.16 37.39 15.03 14.26 12.08 62.21 1.00 1824 1896
## y0[7] 8.08 8.03 15.07 15.30 -16.37 32.22 1.00 1213 1683
## y0[8] 32.22 32.37 14.31 13.92 8.73 56.76 1.00 2053 1973
## y0[9] 36.59 37.03 14.69 14.09 12.26 61.00 1.00 1971 1818
## y0[10] 36.44 36.30 14.57 14.42 12.49 60.49 1.00 2011 2012
## y0[11] 18.92 19.04 14.65 14.13 -5.06 43.58 1.00 1675 1844
## y0[12] 42.43 42.38 15.45 15.16 16.47 67.25 1.00 1844 1847
## y0[13] 43.29 42.90 14.66 14.90 20.46 68.24 1.00 1932 2030
## y0[14] 25.59 25.59 14.69 14.24 1.49 49.29 1.00 1843 1905
## y0[15] 11.71 11.94 15.29 15.18 -13.78 36.49 1.00 1648 1879
## y0[16] 18.60 18.91 14.84 14.60 -5.70 42.02 1.00 1925 1683
## y0[17] 11.28 10.94 15.08 14.58 -13.19 36.95 1.00 1745 1839
## y0[18] 31.31 31.33 14.32 14.34 8.31 55.36 1.00 2056 2006
## y0[19] 8.00 7.99 14.93 14.76 -17.14 33.41 1.00 1034 1705
## y0[20] 32.24 32.20 14.77 14.42 8.33 56.04 1.00 1751 1754
## y0[21] 31.83 31.43 14.72 14.32 7.28 55.90 1.00 1953 1995
## y0[22] 28.43 28.39 14.24 13.99 4.87 51.38 1.00 2016 1934
## y0[23] 43.22 43.49 15.15 14.72 17.60 68.09 1.00 1989 1722
## y0[24] 30.75 30.99 14.71 15.19 6.34 54.50 1.00 1900 2015
## y0[25] 43.02 42.69 15.14 15.02 18.45 68.81 1.00 2072 1875
## y0[26] 33.70 33.65 14.75 14.89 10.40 57.29 1.00 1729 1636
## y0[27] 42.30 42.26 15.36 15.41 16.96 67.53 1.00 2081 1927
## y0[28] 40.85 40.86 15.00 14.69 16.20 65.92 1.00 1869 1757
## y0[29] 14.16 14.16 14.66 14.97 -9.57 39.06 1.00 1640 1910
## y0[30] 8.11 8.26 14.91 14.10 -16.39 33.11 1.00 1476 1577
## y0[31] 25.65 25.86 14.56 14.18 0.96 50.46 1.00 1768 1968
## y0[32] 38.73 39.04 15.43 15.24 13.36 64.20 1.00 1892 2040
## y0[33] 41.40 41.18 14.56 14.27 17.39 65.04 1.00 1958 1711
## y0[34] 14.49 14.03 14.80 15.00 -9.14 39.36 1.00 1714 1916
## y0[35] 21.43 21.23 14.43 14.60 -1.78 45.22 1.00 1775 1930
## y0[36] 9.58 9.75 14.87 14.59 -15.31 34.08 1.00 1481 1891
## y0[37] 30.15 30.34 14.56 14.47 6.20 53.72 1.00 2023 1879
## y0[38] 23.49 23.62 14.72 14.60 -0.71 47.52 1.00 1904 2056
## y0[39] 26.73 26.83 14.77 14.81 2.84 51.30 1.00 2122 1928
## y0[40] 7.48 7.29 15.18 14.82 -17.19 32.84 1.00 1708 1853
## y0[41] 9.83 9.60 14.86 15.40 -14.98 33.65 1.00 1725 1970
## y0[42] 39.37 39.32 14.61 14.13 14.33 64.21 1.00 1798 1932
## y0[43] 19.27 19.17 14.86 14.23 -5.61 42.94 1.00 1676 1971
## y0[44] 22.47 22.54 14.75 14.27 -2.39 46.70 1.00 1997 1880
## y0[45] 35.85 36.18 14.60 14.72 11.69 59.78 1.00 2004 2028
## y0[46] 33.96 34.18 14.72 13.98 9.09 57.90 1.00 1846 2013
## y0[47] 34.48 34.50 14.43 14.13 10.44 57.66 1.00 1957 1952
## y0[48] 33.84 33.63 14.97 14.54 9.09 59.01 1.00 1877 1931
## y0[49] 25.36 25.26 14.42 14.42 1.58 49.08 1.00 1792 1765
## y0[50] 34.96 34.82 14.59 14.02 11.31 59.34 1.00 2067 2196
## m1[1] 7.25 7.27 4.15 3.96 0.89 14.02 1.03 221 275
## m1[2] 11.26 11.32 3.54 3.38 5.74 16.96 1.03 225 294
## m1[3] 15.27 15.30 2.97 2.86 10.55 20.02 1.03 236 305
## m1[4] 19.29 19.30 2.49 2.46 15.31 23.21 1.02 268 332
## m1[5] 23.30 23.34 2.15 2.15 19.91 26.67 1.00 363 388
## m1[6] 27.31 27.33 2.02 2.03 24.11 30.56 1.00 678 757
## m1[7] 31.32 31.31 2.14 2.15 27.86 34.79 1.00 1557 1384
## m1[8] 35.34 35.33 2.47 2.48 31.31 39.24 1.00 1557 1185
## m1[9] 39.35 39.32 2.95 2.91 34.40 44.00 1.00 1082 1076
## m1[10] 43.36 43.32 3.51 3.45 37.58 49.07 1.00 708 936
## y1[1] 6.92 7.19 15.27 15.12 -18.30 31.01 1.00 1204 1764
## y1[2] 10.97 10.99 15.28 14.66 -14.42 36.83 1.00 1168 1708
## y1[3] 15.26 15.51 14.70 14.40 -8.92 38.48 1.00 1360 1595
## y1[4] 19.31 19.33 14.61 13.71 -4.58 43.03 1.00 1684 1822
## y1[5] 23.39 23.62 14.21 13.76 -0.28 46.91 1.00 2110 1920
## y1[6] 27.70 27.36 14.54 13.69 3.73 51.91 1.00 2065 1783
## y1[7] 31.84 32.06 14.92 14.88 7.97 56.48 1.00 2099 1857
## y1[8] 35.15 35.27 14.57 14.38 10.76 58.66 1.00 2019 1785
## y1[9] 39.38 39.16 15.30 14.82 14.23 64.87 1.00 1907 1870
## y1[10] 43.46 43.33 14.89 14.76 18.93 67.66 1.00 2029 1988
estimate as independent class
all b0l,b1l do not have a distribution
class c=1~k ex12-1.stan
yc~N(b0c+b1c*x,s)
data{
int N;
int K;
vector[N] x;
vector[N] y;
array[N] int c;
}
parameters{
vector[K] b0;
vector[K] b1;
real<lower=0> s;
}
model{
for(i in 1:N)
y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0[c[i]]+b1[c[i]]*x[i];
y0[i]=normal_rng(m0[i],s);
}
}mdl=cmdstan_model('./ex12-1.stan')
data=list(N=n,K=k,x=x,y=y,c=factor(c))
mle=mdl$optimize(data=data)## Initial log joint probability = -104749
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmpwb3xsj/model-20a93d7d50fa.stan', line 15, column 4 to column 42)
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 99 -87.4726 0.128499 0.386209 1 1 151
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 199 -87.4689 0.00239103 0.0123312 0.3276 1 264
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 201 -87.4689 0.00222475 0.00169226 1 1 266
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -87.47
## b0[1] 9.92
## b0[2] 4.67
## b0[3] 13.60
## b0[4] -9.46
## b0[5] 13.97
## b0[6] 15.62
## b0[7] 7.27
## b0[8] 9.05
## b0[9] 12.15
## b1[1] 2.81
## b1[2] 0.95
## b1[3] 2.14
## b1[4] 2.02
## b1[5] 2.70
## b1[6] 2.13
## b1[7] 1.57
## b1[8] -0.38
## b1[9] 0.91
## s 3.49
## m0[1] 32.04
## m0[2] 39.81
## m0[3] 5.96
## m0[4] 8.35
## m0[5] 57.91
## m0[6] 2.81
## m0[7] 16.76
## m0[8] 24.86
## m0[9] 26.89
## m0[10] 48.45
## m0[11] 3.65
## m0[12] 30.19
## m0[13] 65.86
## m0[14] 36.91
## m0[15] 7.87
## m0[16] 17.51
## m0[17] 14.49
## m0[18] 24.33
## m0[19] 11.25
## m0[20] 43.58
## m0[21] 24.49
## m0[22] 22.90
## m0[23] 30.35
## m0[24] 16.04
## m0[25] 58.09
## m0[26] 44.59
## m0[27] 56.93
## m0[28] 22.38
## m0[29] 21.95
## m0[30] 15.37
## m0[31] 35.83
## m0[32] 61.39
## m0[33] 22.30
## m0[34] 8.67
## m0[35] 19.23
## m0[36] 16.82
## m0[37] 27.29
## m0[38] 13.11
## m0[39] 22.04
## m0[40] -8.87
## m0[41] 18.66
## m0[42] 25.93
## m0[43] 30.54
## m0[44] 19.94
## m0[45] 56.39
## m0[46] 45.92
## m0[47] 52.31
## m0[48] 54.15
## m0[49] 21.26
## m0[50] 3.21
## y0[1] 31.34
## y0[2] 42.24
## y0[3] 4.15
## y0[4] 11.03
## y0[5] 57.74
## y0[6] 2.90
## y0[7] 18.33
## y0[8] 24.98
## y0[9] 24.09
## y0[10] 45.35
## y0[11] 0.99
## y0[12] 33.97
## y0[13] 72.31
## y0[14] 36.21
## y0[15] 10.63
## y0[16] 16.11
## y0[17] 15.55
## y0[18] 20.10
## y0[19] 12.28
## y0[20] 41.61
## y0[21] 24.10
## y0[22] 21.15
## y0[23] 32.76
## y0[24] 18.54
## y0[25] 58.97
## y0[26] 38.41
## y0[27] 59.37
## y0[28] 20.33
## y0[29] 22.44
## y0[30] 19.07
## y0[31] 32.60
## y0[32] 59.67
## y0[33] 24.00
## y0[34] 11.27
## y0[35] 17.70
## y0[36] 21.32
## y0[37] 27.99
## y0[38] 8.81
## y0[39] 25.00
## y0[40] -13.42
## y0[41] 21.28
## y0[42] 30.65
## y0[43] 32.41
## y0[44] 15.24
## y0[45] 55.36
## y0[46] 50.97
## y0[47] 54.41
## y0[48] 55.49
## y0[49] 18.54
## y0[50] 0.42mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 1.3 seconds.
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 finished in 1.5 seconds.
## Chain 4 finished in 1.4 seconds.
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 1.7 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 1.5 seconds.
## Total execution time: 1.8 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -98.10 -97.65 3.89 3.74 -105.02 -92.64 1.00 623 1086
## b0[1] 9.93 9.87 4.62 4.70 2.54 17.27 1.00 1925 1389
## b0[2] 4.76 4.75 5.00 4.85 -3.24 12.98 1.00 1602 1272
## b0[3] 13.63 13.68 2.99 2.94 8.87 18.54 1.00 2880 1570
## b0[4] -9.48 -9.39 3.80 3.68 -15.76 -3.18 1.00 2324 1476
## b0[5] 14.13 14.11 4.90 4.98 5.92 22.38 1.00 1801 1257
## b0[6] 15.62 15.65 3.73 3.50 9.40 21.79 1.00 2324 1364
## b0[7] 7.30 7.31 4.12 4.16 0.66 13.93 1.00 1923 1104
## b0[8] 9.02 9.09 5.24 5.37 0.22 17.24 1.00 1407 1275
## b0[9] 12.19 12.08 4.29 4.36 5.44 19.30 1.00 1951 1363
## b1[1] 2.82 2.82 0.32 0.31 2.29 3.34 1.00 1946 1221
## b1[2] 0.94 0.95 0.35 0.35 0.35 1.51 1.00 1677 1597
## b1[3] 2.14 2.13 0.26 0.26 1.71 2.57 1.00 2581 1653
## b1[4] 2.02 2.02 0.26 0.24 1.60 2.44 1.00 2237 1306
## b1[5] 2.68 2.68 0.36 0.36 2.10 3.26 1.00 1850 1398
## b1[6] 2.13 2.12 0.26 0.24 1.72 2.55 1.00 2389 1685
## b1[7] 1.57 1.57 0.39 0.39 0.93 2.19 1.00 1951 1477
## b1[8] -0.38 -0.38 0.43 0.42 -1.05 0.35 1.01 1566 1441
## b1[9] 0.90 0.91 0.38 0.38 0.28 1.52 1.00 1949 1617
## s 4.53 4.46 0.60 0.57 3.65 5.63 1.00 1078 1349
## m0[1] 32.10 32.10 3.42 3.49 26.59 37.69 1.00 1910 1704
## m0[2] 39.83 39.84 1.61 1.58 37.21 42.43 1.00 2191 1232
## m0[3] 5.93 6.06 2.59 2.58 1.65 10.03 1.00 1389 1490
## m0[4] 8.38 8.40 3.90 3.93 2.11 14.70 1.00 1916 1102
## m0[5] 57.86 57.79 2.65 2.52 53.56 62.29 1.00 2124 1283
## m0[6] 2.79 2.71 3.24 3.17 -2.36 8.17 1.00 1958 1365
## m0[7] 16.90 16.92 4.58 4.66 9.32 24.62 1.00 1804 1258
## m0[8] 24.84 24.86 1.92 1.79 21.75 28.00 1.00 2028 1611
## m0[9] 26.85 26.83 2.57 2.43 22.66 30.94 1.00 1996 1585
## m0[10] 48.47 48.49 2.26 2.24 44.76 52.09 1.00 2195 1357
## m0[11] 3.62 3.66 2.47 2.41 -0.50 7.58 1.00 2256 1364
## m0[12] 30.15 30.17 2.55 2.44 26.10 34.40 1.00 1954 1363
## m0[13] 66.03 66.03 3.68 3.68 60.00 72.00 1.01 2038 1522
## m0[14] 36.88 36.94 2.01 1.94 33.54 40.11 1.00 2301 1512
## m0[15] 7.84 7.85 4.08 4.12 0.89 14.42 1.00 1388 1281
## m0[16] 17.56 17.56 2.41 2.35 13.58 21.50 1.00 1835 1202
## m0[17] 14.51 14.39 3.38 3.38 9.20 20.03 1.00 1938 1283
## m0[18] 24.30 24.30 1.78 1.70 21.47 27.27 1.00 2044 1560
## m0[19] 11.27 11.22 4.49 4.60 4.06 18.41 1.00 1929 1410
## m0[20] 43.60 43.63 1.85 1.77 40.57 46.58 1.00 2144 1218
## m0[21] 24.47 24.45 1.82 1.74 21.53 27.46 1.00 2038 1613
## m0[22] 22.88 22.85 1.50 1.45 20.39 25.31 1.00 2104 1571
## m0[23] 30.31 30.31 2.57 2.45 26.24 34.58 1.00 1953 1363
## m0[24] 16.00 16.02 1.84 1.78 13.01 19.00 1.00 2011 1460
## m0[25] 58.04 57.96 2.66 2.54 53.72 62.49 1.00 2125 1283
## m0[26] 44.61 44.64 1.93 1.87 41.44 47.70 1.00 2149 1233
## m0[27] 56.89 56.84 2.57 2.42 52.69 61.21 1.00 2120 1251
## m0[28] 22.37 22.32 3.09 3.27 17.45 27.46 1.00 2030 1751
## m0[29] 21.98 21.96 2.16 2.17 18.58 25.54 1.00 2891 1621
## m0[30] 15.40 15.44 2.81 2.74 10.93 20.00 1.00 2890 1523
## m0[31] 35.86 35.85 1.49 1.41 33.41 38.28 1.00 2296 1364
## m0[32] 61.25 61.14 3.04 2.95 56.30 66.31 1.00 1969 1409
## m0[33] 22.29 22.25 3.07 3.24 17.40 27.33 1.00 2032 1751
## m0[34] 8.74 8.76 3.73 3.52 2.74 14.82 1.00 1612 1365
## m0[35] 19.24 19.19 1.79 1.74 16.40 22.19 1.00 1941 1619
## m0[36] 16.85 16.89 2.66 2.61 12.55 21.18 1.00 2890 1522
## m0[37] 27.35 27.34 2.64 2.65 23.13 31.64 1.00 1872 1414
## m0[38] 13.15 13.16 2.60 2.46 8.69 17.45 1.00 1721 1605
## m0[39] 22.03 22.02 1.43 1.41 19.71 24.35 1.01 2101 1542
## m0[40] -8.90 -8.82 3.74 3.62 -15.08 -2.72 1.00 2322 1406
## m0[41] 18.65 18.63 3.42 3.23 12.98 24.29 1.00 2315 1325
## m0[42] 25.88 25.88 2.21 2.13 22.30 29.56 1.00 1925 1487
## m0[43] 30.52 30.56 2.37 2.24 26.59 34.45 1.00 2282 1469
## m0[44] 19.94 19.90 1.62 1.64 17.38 22.61 1.00 1963 1682
## m0[45] 56.29 56.28 2.67 2.56 51.98 60.68 1.00 1988 1416
## m0[46] 45.94 45.97 2.03 1.97 42.57 49.17 1.00 2166 1326
## m0[47] 52.45 52.48 2.79 2.73 47.95 56.92 1.01 2051 1521
## m0[48] 54.06 54.06 2.54 2.44 50.03 58.26 1.00 1998 1375
## m0[49] 21.26 21.26 1.44 1.45 18.94 23.59 1.00 2046 1662
## m0[50] 3.19 3.14 2.94 2.90 -1.52 8.07 1.00 1981 1464
## y0[1] 32.05 32.05 5.71 5.60 22.58 41.35 1.00 1840 1920
## y0[2] 39.87 39.96 4.74 4.72 31.90 47.42 1.00 2021 1618
## y0[3] 5.92 5.90 5.24 5.15 -2.66 14.73 1.00 1793 1895
## y0[4] 8.43 8.43 5.98 5.93 -1.30 18.13 1.00 1993 1710
## y0[5] 57.86 57.82 5.47 5.26 49.06 67.06 1.00 1925 1744
## y0[6] 2.83 2.74 5.65 5.61 -6.39 12.38 1.00 1974 1760
## y0[7] 16.87 16.93 6.46 6.26 5.85 27.24 1.00 2042 1851
## y0[8] 24.77 24.73 4.88 4.85 16.74 32.92 1.00 1517 1820
## y0[9] 26.83 26.82 5.23 5.20 18.32 35.32 1.00 2029 1899
## y0[10] 48.50 48.53 4.91 4.85 40.46 56.59 1.00 2025 2013
## y0[11] 3.68 3.81 5.20 5.14 -5.10 12.10 1.00 1937 1729
## y0[12] 30.17 30.10 5.28 5.25 21.39 38.97 1.00 1864 1720
## y0[13] 65.93 65.83 5.76 5.65 56.41 75.70 1.00 1869 1639
## y0[14] 36.86 36.85 5.08 5.18 28.63 45.18 1.00 2238 2099
## y0[15] 7.81 7.98 5.94 5.76 -1.97 17.25 1.00 1725 1622
## y0[16] 17.58 17.65 5.10 4.98 9.13 25.83 1.00 1850 1782
## y0[17] 14.52 14.57 5.63 5.72 5.43 23.97 1.00 1866 1702
## y0[18] 24.42 24.35 4.73 4.71 16.67 32.08 1.00 1488 1684
## y0[19] 11.37 11.45 6.17 5.82 1.37 21.70 1.00 1973 1701
## y0[20] 43.64 43.52 4.86 4.70 35.95 51.34 1.00 2087 1970
## y0[21] 24.44 24.37 4.85 4.77 16.54 32.33 1.00 1843 1891
## y0[22] 22.86 22.92 4.71 4.66 14.95 30.46 1.00 1917 1971
## y0[23] 30.25 30.20 5.21 5.09 22.02 38.80 1.00 1965 1868
## y0[24] 15.99 15.81 5.03 4.79 8.25 24.69 1.00 1650 1912
## y0[25] 58.07 58.08 5.40 5.38 49.07 67.23 1.00 2117 1539
## y0[26] 44.69 44.59 4.96 5.01 36.62 52.94 1.00 2033 1800
## y0[27] 56.73 56.75 5.33 5.09 47.97 65.14 1.00 1994 1945
## y0[28] 22.35 22.25 5.56 5.29 13.55 31.40 1.00 2009 1911
## y0[29] 22.02 22.03 5.06 5.06 14.10 30.20 1.00 2136 1808
## y0[30] 15.41 15.41 5.40 5.16 6.35 24.08 1.00 2181 1880
## y0[31] 35.81 35.84 4.85 4.92 27.83 43.79 1.00 1883 1872
## y0[32] 61.12 61.21 5.40 5.34 52.21 69.92 1.00 1838 1606
## y0[33] 22.41 22.57 5.58 5.57 13.08 31.42 1.00 2058 1848
## y0[34] 8.72 8.55 6.05 6.05 -0.99 18.98 1.00 1811 1729
## y0[35] 19.06 19.08 4.90 4.82 10.98 27.15 1.00 2170 1895
## y0[36] 16.90 16.82 5.26 5.11 7.84 25.43 1.00 2253 1879
## y0[37] 27.20 27.18 5.22 5.04 18.80 36.18 1.00 2081 1921
## y0[38] 13.15 13.17 5.23 5.17 4.54 21.73 1.00 2053 1754
## y0[39] 22.03 22.10 4.79 4.64 14.35 29.89 1.00 1873 1949
## y0[40] -8.95 -9.00 5.86 5.72 -18.40 0.48 1.00 2182 1801
## y0[41] 18.67 18.56 5.73 5.63 9.31 28.24 1.00 2071 1624
## y0[42] 25.94 25.83 5.03 4.96 17.72 34.26 1.00 2074 1835
## y0[43] 30.53 30.57 5.10 5.02 22.31 39.07 1.00 1991 1877
## y0[44] 19.98 19.96 4.86 4.83 12.03 27.89 1.00 2031 1767
## y0[45] 56.36 56.35 5.30 5.04 47.58 65.06 1.00 1862 1808
## y0[46] 45.75 45.70 4.97 4.72 37.83 53.93 1.00 2176 2054
## y0[47] 52.39 52.33 5.38 5.34 43.84 61.22 1.00 1978 1932
## y0[48] 54.02 54.01 5.13 5.08 45.56 62.22 1.00 1993 1892
## y0[49] 21.33 21.28 4.64 4.38 13.63 29.16 1.00 1920 1868
## y0[50] 3.25 3.27 5.41 5.28 -5.61 12.52 1.00 2043 2099
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,shape=c,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
#geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
#geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
#geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
#geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
qplot(x,y,facets=~c,size=I(2),col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
estimate as class have relation
all b0l,b1l have a distribution
class c=1~k
b0c~N(b00,sb0)
b1c~N(b10,sb1)
yc~N(b0c+b1c*x,s)ex12-2.stan
class have relation
data{
int N;
int K;
vector[N] x;
vector[N] y;
array[N] int c;
}
parameters{
real b00;
real<lower=0,upper=100> sb0;
vector[K] b0;
real b10;
real<lower=0,upper=100> sb1;
vector[K] b1;
real<lower=0,upper=100> s;
}
model{
b0~normal(b00,sb0);
b1~normal(b10,sb1);
for(i in 1:N)
y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0[c[i]]+b1[c[i]]*x[i];
y0[i]=normal_rng(m0[i],s);
}
}mdl=cmdstan_model('./ex12-2.stan')
data=list(N=n,K=k,x=x,y=y,c=factor(c))
mle=mdl$optimize(data=data)## Initial log joint probability = -314.534
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 99 -36.5232 0.301266 99672.6 0.141 0.77 139
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 199 121.716 0.391023 1.47123e+14 0.08681 0.5203 311
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 220 139.916 0.289332 4.96632e+15 1e-12 0.001 406 LS failed, Hessian reset
## Finished in 0.1 seconds.try(print(mle))## Error : Fitting failed. Unable to print.mcmc=goMCMC(mdl,data,wrm=1000,smp=2000)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 3 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 1 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 1 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 2 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 2 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 3 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 3 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 4 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 4 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 1 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 2 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 3 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 4 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 1 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 2 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 4 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 1 finished in 1.3 seconds.
## Chain 2 finished in 1.3 seconds.
## Chain 4 finished in 1.3 seconds.
## Chain 3 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 3 finished in 1.3 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 1.3 seconds.
## Total execution time: 1.5 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -123.72 -123.26 4.24 4.07 -131.31 -117.68 1.00 2474 4034
## b00 8.45 8.46 3.40 3.04 2.99 13.78 1.00 7156 5013
## sb0 8.62 8.02 3.35 2.76 4.55 14.78 1.00 4184 4309
## b0[1] 10.44 10.43 3.96 3.88 3.99 17.11 1.00 8744 6284
## b0[2] 5.03 5.03 4.31 4.22 -1.98 12.14 1.00 7276 5822
## b0[3] 13.11 13.15 2.80 2.74 8.55 17.71 1.00 6169 5776
## b0[4] -5.80 -5.88 3.99 3.95 -12.22 0.86 1.00 6285 4906
## b0[5] 13.44 13.35 4.13 4.10 6.84 20.35 1.00 7094 6068
## b0[6] 14.48 14.47 3.51 3.47 8.73 20.25 1.00 6966 5527
## b0[7] 7.40 7.42 3.61 3.47 1.47 13.38 1.00 8428 5679
## b0[8] 6.74 6.83 4.32 4.17 -0.60 13.67 1.00 7687 6211
## b0[9] 10.50 10.54 3.77 3.65 4.24 16.64 1.00 8060 6066
## b10 1.67 1.68 0.43 0.38 0.98 2.36 1.00 6651 4829
## sb1 1.18 1.10 0.43 0.32 0.68 1.97 1.00 6527 4576
## b1[1] 2.74 2.75 0.29 0.29 2.26 3.22 1.00 8492 6055
## b1[2] 0.95 0.94 0.32 0.31 0.43 1.47 1.00 7462 6178
## b1[3] 2.17 2.17 0.25 0.24 1.77 2.58 1.00 7221 5744
## b1[4] 1.80 1.81 0.27 0.26 1.35 2.23 1.00 6671 5405
## b1[5] 2.71 2.72 0.30 0.30 2.20 3.21 1.00 7494 6070
## b1[6] 2.20 2.20 0.24 0.23 1.81 2.59 1.00 7460 5214
## b1[7] 1.56 1.56 0.34 0.34 1.00 2.12 1.00 8377 5988
## b1[8] -0.14 -0.15 0.37 0.36 -0.74 0.49 1.00 7696 5736
## b1[9] 1.05 1.05 0.33 0.33 0.52 1.60 1.00 8216 6128
## s 4.52 4.45 0.59 0.57 3.69 5.60 1.00 5091 5985
## m0[1] 32.02 32.00 3.40 3.35 26.38 37.60 1.00 8126 6238
## m0[2] 39.68 39.70 1.63 1.63 37.04 42.33 1.00 8875 5062
## m0[3] 5.57 5.59 2.30 2.23 1.78 9.30 1.00 8821 6448
## m0[4] 8.47 8.49 3.43 3.30 2.80 14.16 1.00 8413 5868
## m0[5] 58.02 57.99 2.52 2.54 53.94 62.16 1.00 8305 6408
## m0[6] 4.39 4.33 3.27 3.24 -0.93 9.87 1.00 8504 6460
## m0[7] 16.23 16.14 3.87 3.82 10.05 22.72 1.00 7094 6361
## m0[8] 25.25 25.27 1.84 1.80 22.24 28.28 1.00 8373 6156
## m0[9] 27.60 27.59 2.37 2.32 23.78 31.50 1.00 8430 6514
## m0[10] 48.44 48.44 2.24 2.23 44.79 52.14 1.00 8447 6094
## m0[11] 5.88 5.83 2.60 2.54 1.71 10.17 1.00 6489 5230
## m0[12] 29.56 29.59 2.57 2.51 25.30 33.79 1.00 8838 6296
## m0[13] 65.12 65.18 3.72 3.65 58.97 71.22 1.00 8827 6172
## m0[14] 36.40 36.43 1.93 1.91 33.23 39.56 1.00 7800 5755
## m0[15] 6.29 6.37 3.37 3.25 0.58 11.72 1.00 7783 6676
## m0[16] 17.58 17.57 2.30 2.27 13.81 21.36 1.00 8101 6628
## m0[17] 13.21 13.23 3.00 2.90 8.28 18.10 1.00 8043 5919
## m0[18] 24.62 24.65 1.73 1.69 21.80 27.47 1.00 8338 6055
## m0[19] 11.75 11.72 3.86 3.78 5.45 18.23 1.00 8760 6389
## m0[20] 43.50 43.51 1.86 1.83 40.48 46.54 1.00 8703 5623
## m0[21] 24.81 24.84 1.76 1.72 21.93 27.70 1.00 8349 6228
## m0[22] 22.97 22.96 1.51 1.46 20.49 25.40 1.00 8237 5028
## m0[23] 29.70 29.73 2.59 2.53 25.42 33.94 1.00 8837 6310
## m0[24] 16.94 16.90 1.91 1.86 13.82 20.10 1.00 7975 5092
## m0[25] 58.21 58.18 2.54 2.55 54.10 62.37 1.00 8302 6173
## m0[26] 44.52 44.53 1.93 1.92 41.39 47.67 1.00 8647 5684
## m0[27] 57.02 57.01 2.45 2.45 53.04 61.02 1.00 8319 6214
## m0[28] 22.71 22.65 3.08 3.00 17.78 27.74 1.00 8412 6197
## m0[29] 21.57 21.58 2.04 2.02 18.23 24.96 1.00 6358 5809
## m0[30] 14.90 14.94 2.62 2.58 10.63 19.22 1.00 6173 5631
## m0[31] 35.65 35.67 1.51 1.50 33.20 38.10 1.00 8852 5134
## m0[32] 61.02 61.01 2.86 2.83 56.36 65.69 1.00 8262 5445
## m0[33] 22.63 22.57 3.06 2.99 17.73 27.64 1.00 8420 6097
## m0[34] 9.02 9.04 3.23 3.17 3.78 14.40 1.00 7418 6338
## m0[35] 18.72 18.74 1.70 1.65 15.91 21.46 1.00 7850 6647
## m0[36] 16.37 16.40 2.49 2.43 12.35 20.45 1.00 6193 5630
## m0[37] 27.31 27.28 2.70 2.67 22.88 31.74 1.00 8025 5571
## m0[38] 13.45 13.46 2.36 2.31 9.62 17.35 1.00 7962 6548
## m0[39] 21.97 21.97 1.45 1.39 19.57 24.28 1.00 8209 5373
## m0[40] -5.29 -5.36 3.92 3.88 -11.57 1.27 1.00 6289 4942
## m0[41] 17.61 17.58 3.22 3.19 12.30 22.91 1.00 6942 5728
## m0[42] 25.75 25.76 2.23 2.16 22.05 29.45 1.00 8854 5711
## m0[43] 29.84 29.85 2.26 2.22 26.15 33.54 1.00 7298 5730
## m0[44] 19.54 19.57 1.58 1.53 16.93 22.10 1.00 7910 6227
## m0[45] 56.01 55.99 2.54 2.49 51.86 60.13 1.00 8286 5270
## m0[46] 45.87 45.87 2.03 2.01 42.58 49.22 1.00 8577 5753
## m0[47] 51.88 51.92 2.86 2.81 47.14 56.52 1.00 9043 5148
## m0[48] 53.75 53.74 2.42 2.37 49.81 57.69 1.00 8283 5458
## m0[49] 21.07 21.09 1.45 1.41 18.68 23.40 1.00 8210 5651
## m0[50] 4.54 4.51 2.99 2.92 -0.37 9.52 1.00 8633 6172
## y0[1] 31.97 31.91 5.66 5.55 22.79 41.23 1.00 7778 7766
## y0[2] 39.77 39.76 4.87 4.76 31.73 47.89 1.00 8342 7861
## y0[3] 5.61 5.57 5.07 4.92 -2.66 13.89 1.00 8037 7613
## y0[4] 8.50 8.54 5.66 5.60 -0.79 17.62 1.00 7871 7393
## y0[5] 58.03 58.02 5.18 5.10 49.55 66.53 1.00 8243 7442
## y0[6] 4.39 4.31 5.66 5.57 -4.62 13.75 1.00 7883 7402
## y0[7] 16.17 16.24 6.06 5.78 6.27 26.27 1.00 6878 7230
## y0[8] 25.19 25.29 4.98 4.80 16.97 33.20 1.00 8209 7187
## y0[9] 27.49 27.47 5.14 5.06 18.97 35.99 1.00 7292 7271
## y0[10] 48.47 48.48 5.09 4.89 40.12 56.94 1.00 8290 7837
## y0[11] 5.85 5.72 5.24 5.18 -2.63 14.52 1.00 7830 7396
## y0[12] 29.52 29.45 5.25 5.09 20.91 38.23 1.00 8493 7789
## y0[13] 65.17 65.26 5.90 5.73 55.52 74.73 1.00 8230 7849
## y0[14] 36.35 36.40 4.99 4.85 28.14 44.55 1.00 7970 7757
## y0[15] 6.34 6.37 5.64 5.54 -3.11 15.46 1.00 8018 7295
## y0[16] 17.60 17.51 5.11 5.01 9.39 25.99 1.00 7777 7263
## y0[17] 13.20 13.19 5.40 5.34 4.33 22.12 1.00 7789 7808
## y0[18] 24.69 24.69 4.92 4.86 16.68 32.83 1.00 7981 8098
## y0[19] 11.78 11.85 6.09 6.05 1.52 21.65 1.00 8346 7375
## y0[20] 43.55 43.51 4.93 4.86 35.53 51.58 1.00 8138 7794
## y0[21] 24.89 24.94 4.89 4.82 16.82 32.78 1.00 8370 7850
## y0[22] 23.03 23.03 4.79 4.69 15.26 30.87 1.00 7878 7842
## y0[23] 29.77 29.77 5.22 5.14 21.27 38.39 1.00 8240 8098
## y0[24] 16.89 16.83 4.96 4.91 8.78 24.99 1.00 8019 7504
## y0[25] 58.24 58.27 5.23 5.12 49.66 66.84 1.00 8108 7101
## y0[26] 44.67 44.68 4.82 4.71 36.84 52.53 1.00 7954 7626
## y0[27] 57.02 57.05 5.17 5.12 48.60 65.50 1.00 7939 7647
## y0[28] 22.72 22.72 5.50 5.52 13.70 31.73 1.00 7759 7674
## y0[29] 21.59 21.69 5.03 4.92 13.31 29.63 1.00 7534 7540
## y0[30] 14.93 15.04 5.25 5.17 6.12 23.40 1.00 7651 7808
## y0[31] 35.68 35.66 4.79 4.75 27.82 43.48 1.00 8558 7674
## y0[32] 61.03 60.99 5.34 5.21 52.29 69.81 1.00 7870 7778
## y0[33] 22.67 22.65 5.47 5.41 13.83 31.72 1.00 8005 7076
## y0[34] 8.96 8.86 5.55 5.50 -0.10 18.10 1.00 7593 7594
## y0[35] 18.69 18.73 4.92 4.74 10.53 26.70 1.00 7900 7498
## y0[36] 16.35 16.36 5.23 5.16 7.65 24.92 1.00 7336 7398
## y0[37] 27.31 27.26 5.32 5.23 18.60 36.15 1.00 7784 7633
## y0[38] 13.43 13.39 5.14 5.04 5.08 22.02 1.00 8160 7544
## y0[39] 21.97 21.97 4.77 4.62 14.13 29.86 1.00 7064 7320
## y0[40] -5.23 -5.26 6.06 5.91 -15.09 4.83 1.00 7198 6895
## y0[41] 17.55 17.57 5.61 5.48 8.39 26.85 1.00 8214 7125
## y0[42] 25.72 25.67 5.10 5.03 17.53 34.16 1.00 8198 7550
## y0[43] 29.85 29.83 5.10 5.06 21.52 38.16 1.00 7791 7972
## y0[44] 19.56 19.60 4.80 4.67 11.62 27.42 1.00 8096 7875
## y0[45] 55.97 56.09 5.21 5.11 47.45 64.32 1.00 8078 7710
## y0[46] 45.84 45.77 4.94 4.88 37.78 54.10 1.00 7572 7481
## y0[47] 51.85 51.88 5.41 5.22 43.04 60.84 1.00 8736 7203
## y0[48] 53.68 53.74 5.13 5.03 45.32 62.07 1.00 8136 6479
## y0[49] 21.05 21.00 4.76 4.66 13.27 28.92 1.00 7153 7071
## y0[50] 4.62 4.49 5.47 5.37 -4.18 13.72 1.00 8256 7782
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,shape=c,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
#geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
#geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
#geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
#geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
qplot(x,y,facets=~c,size=I(2),col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
generalized linear mixed model
(X,y)i=1-n
b[b0,b1,...]
linear model y~N(Xb,s)
generalized linear model y~dist.(m=link(Xb),s)
fixed effect b0, b1,...
individual random effect b0+r0i r0~N(0,sr0), b1+r1i r1~N(0,sr1),...
class c=1-k
class effect b0+r0c r0~N(0,sr0), b1+r1c r1~N(0,sr1),...
for y=dist.(m,s)
random intercept model m=b0+r0i+b1*x, m=b0+r0c+b1*x
random coefficient model m=b0+(b1+r1i)*x, m=b0+(b1+r1c)*x
mixed model m=b0+r0i+(b1+r1i)*x, m=b0+r0c+(b1+r1c)*x
note
@ yi=b0+b1*xi+r0i is not useful, ri is included in s
@ yi=b0+b1*xi+r0c is useful, repeated observation can be treated by class effect
@ when appling Poisson dist.(y is larger than 0, error is larger at large x),
but variance is larger than mean (over dispersion),
random intercept model is usefuln=20
x=runif(n,0,10)
r0=rnorm(n,0,1)
y=rpois(n,exp(2+0.1*x+r0))
qplot(x,y)
ex13-0.stan
generalized linear model, poisson regression
data{
int N;
vector[N] x;
array[N] int y;
}
parameters{
real b0;
real b1;
}
model{
y~poisson_log(b0+b1*x);
}
generated quantities{
vector[N] m0;
array[N] int y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=poisson_log_rng(m0[i]);
}
}mdl=cmdstan_model('./ex13-0.stan')
data=list(N=n,x=x,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -8.45564e+06
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 21 1184.69 0.000250358 0.0562033 0.9629 0.9629 61
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 1184.69
## b0 2.80
## b1 0.08
## m0[1] 3.19
## m0[2] 3.26
## m0[3] 3.56
## m0[4] 3.32
## m0[5] 3.31
## m0[6] 3.43
## m0[7] 3.30
## m0[8] 3.11
## m0[9] 3.40
## m0[10] 2.82
## m0[11] 2.81
## m0[12] 3.46
## m0[13] 3.46
## m0[14] 2.98
## m0[15] 2.86
## m0[16] 3.28
## m0[17] 3.15
## m0[18] 3.52
## m0[19] 3.50
## m0[20] 2.86
## y0[1] 19.00
## y0[2] 23.00
## y0[3] 31.00
## y0[4] 24.00
## y0[5] 34.00
## y0[6] 30.00
## y0[7] 32.00
## y0[8] 20.00
## y0[9] 38.00
## y0[10] 13.00
## y0[11] 17.00
## y0[12] 29.00
## y0[13] 33.00
## y0[14] 21.00
## y0[15] 24.00
## y0[16] 21.00
## y0[17] 32.00
## y0[18] 34.00
## y0[19] 34.00
## y0[20] 17.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 1183.21 1183.95 4.12 0.77 1181.28 1184.64 1.01 243 148
## b0 2.78 2.79 0.15 0.11 2.60 2.96 1.01 276 122
## b1 0.08 0.08 0.02 0.02 0.05 0.11 1.01 297 138
## m0[1] 3.18 3.19 0.06 0.05 3.10 3.26 1.01 374 147
## m0[2] 3.25 3.25 0.05 0.05 3.17 3.33 1.01 461 159
## m0[3] 3.56 3.56 0.08 0.07 3.43 3.68 1.01 495 266
## m0[4] 3.32 3.32 0.05 0.05 3.24 3.39 1.01 671 243
## m0[5] 3.30 3.30 0.05 0.05 3.22 3.38 1.01 594 238
## m0[6] 3.42 3.43 0.05 0.05 3.33 3.51 1.00 1643 1532
## m0[7] 3.30 3.30 0.05 0.05 3.22 3.37 1.01 570 236
## m0[8] 3.11 3.11 0.08 0.06 3.01 3.20 1.01 319 127
## m0[9] 3.40 3.40 0.05 0.05 3.31 3.48 1.00 1873 1457
## m0[10] 2.80 2.81 0.15 0.11 2.62 2.97 1.01 276 120
## m0[11] 2.79 2.81 0.15 0.11 2.61 2.97 1.01 276 120
## m0[12] 3.46 3.46 0.06 0.06 3.36 3.56 1.00 982 666
## m0[13] 3.46 3.46 0.06 0.06 3.36 3.56 1.00 974 666
## m0[14] 2.97 2.98 0.11 0.08 2.85 3.10 1.01 286 123
## m0[15] 2.84 2.85 0.14 0.10 2.67 3.00 1.01 277 123
## m0[16] 3.27 3.27 0.05 0.05 3.19 3.35 1.01 503 179
## m0[17] 3.14 3.15 0.07 0.05 3.05 3.23 1.01 340 116
## m0[18] 3.52 3.52 0.07 0.06 3.40 3.62 1.01 580 322
## m0[19] 3.50 3.50 0.07 0.06 3.38 3.60 1.01 655 391
## m0[20] 2.85 2.86 0.14 0.10 2.68 3.01 1.01 277 121
## y0[1] 24.02 24.00 5.21 5.93 16.00 33.00 1.00 1493 1602
## y0[2] 25.83 26.00 5.24 5.93 17.00 35.00 1.00 1743 1864
## y0[3] 35.38 35.00 6.50 5.93 25.00 47.00 1.00 1541 1813
## y0[4] 27.68 27.00 5.30 5.93 19.00 36.05 1.00 1906 2038
## y0[5] 27.13 27.00 5.60 5.93 18.00 36.00 1.00 1691 1828
## y0[6] 30.72 30.50 5.71 5.19 22.00 41.00 1.00 1968 1722
## y0[7] 27.18 27.00 5.36 5.93 19.00 36.00 1.00 1876 1885
## y0[8] 22.36 22.00 5.05 4.45 14.00 31.00 1.01 1015 1039
## y0[9] 29.89 30.00 5.64 5.93 21.00 40.00 1.00 1969 1858
## y0[10] 16.59 16.00 4.61 4.45 9.00 24.00 1.00 650 463
## y0[11] 16.39 16.00 4.54 4.45 9.00 24.00 1.00 630 358
## y0[12] 31.80 32.00 6.02 5.93 22.00 42.00 1.00 2125 1941
## y0[13] 31.87 32.00 5.79 5.93 23.00 42.00 1.00 1996 1941
## y0[14] 19.52 19.00 4.83 4.45 12.00 28.00 1.00 846 511
## y0[15] 17.25 17.00 4.71 4.45 10.00 25.00 1.00 654 493
## y0[16] 26.39 26.00 5.38 5.93 18.00 36.00 1.00 1722 1847
## y0[17] 23.11 23.00 5.04 4.45 15.00 31.00 1.00 1185 923
## y0[18] 33.70 34.00 6.26 5.93 24.00 44.00 1.00 1955 1254
## y0[19] 32.93 33.00 6.06 5.93 23.00 43.00 1.00 1276 1800
## y0[20] 17.19 17.00 4.52 4.45 10.00 25.00 1.00 650 401
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=exp(m)),xy,col='black')
ex13-1.stan
generalized linear mixed model, poisson regression
data{
int N;
vector[N] x;
array[N] int y;
}
parameters{
real b0;
real b1;
real<lower=0> sr0;
vector[N] r0;
}
model{
r0~normal(0,sr0);
for(i in 1:N)
y~poisson_log(b0+r0[i]+b1*x);
}
generated quantities{
vector[N] m0;
array[N] int y0;
for(i in 1:N){
m0[i]=b0+r0[i]+b1*x[i];
y0[i]=poisson_log_rng(m0[i]);
}
}mdl=cmdstan_model('./ex13-1.stan')
data=list(N=n,x=x,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -7.77186e+07
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 99 23721.8 0.404198 109.714 0.9746 0.9746 153
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 199 24031.5 0.425447 1.22229e+08 0.1516 1 314
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 299 24418.7 0.80365 9.7267e+18 1.741 0.5183 481
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 399 24678.1 0.0183169 1.03074e+29 0.09589 0.09589 669
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 499 24817.7 0.2255 3.83185e+39 2.326 0.3774 857
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 541 24829.3 0.044588 2.15164e+41 0.8437 0.8437 923
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 24829.30
## b0 2.51
## b1 0.17
## sr0 0.00
## r0[1] 0.00
## r0[2] 0.00
## r0[3] 0.00
## r0[4] 0.00
## r0[5] 0.00
## r0[6] 0.00
## r0[7] 0.00
## r0[8] 0.00
## r0[9] 0.00
## r0[10] 0.00
## r0[11] 0.00
## r0[12] 0.00
## r0[13] 0.00
## r0[14] 0.00
## r0[15] 0.00
## r0[16] 0.00
## r0[17] 0.00
## r0[18] 0.00
## r0[19] 0.00
## r0[20] 0.00
## m0[1] 3.36
## m0[2] 3.51
## m0[3] 4.17
## m0[4] 3.65
## m0[5] 3.62
## m0[6] 3.87
## m0[7] 3.61
## m0[8] 3.20
## m0[9] 3.81
## m0[10] 2.56
## m0[11] 2.54
## m0[12] 3.95
## m0[13] 3.95
## m0[14] 2.92
## m0[15] 2.64
## m0[16] 3.55
## m0[17] 3.28
## m0[18] 4.07
## m0[19] 4.02
## m0[20] 2.65
## y0[1] 24.00
## y0[2] 32.00
## y0[3] 63.00
## y0[4] 34.00
## y0[5] 41.00
## y0[6] 41.00
## y0[7] 40.00
## y0[8] 23.00
## y0[9] 58.00
## y0[10] 14.00
## y0[11] 13.00
## y0[12] 50.00
## y0[13] 52.00
## y0[14] 17.00
## y0[15] 8.00
## y0[16] 30.00
## y0[17] 30.00
## y0[18] 68.00
## y0[19] 61.00
## y0[20] 9.00mcmc=goMCMC(mdl,data,wrm=500)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 3 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 3 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 4 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 1 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 2 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 1 finished in 1.2 seconds.
## Chain 2 finished in 1.1 seconds.
## Chain 3 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 3 finished in 1.3 seconds.
## Chain 4 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 4 finished in 1.3 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 1.2 seconds.
## Total execution time: 1.5 seconds.seeMCMC(mcmc,'[m,r]')## variable mean median sd mad q5 q95 rhat ess_bulk
## lp__ 23773.49 23771.70 14.54 14.23 23752.20 23799.80 1.08 51
## b0 2.80 2.79 0.02 0.02 2.75 2.83 1.01 768
## b1 0.08 0.08 0.00 0.00 0.07 0.08 1.01 836
## sr0 0.01 0.01 0.01 0.01 0.00 0.02 1.07 59
## r0[1] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 3212
## r0[2] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 2708
## r0[3] 0.00 0.00 0.01 0.01 -0.02 0.02 1.04 2613
## r0[4] 0.00 0.00 0.01 0.01 -0.02 0.02 1.01 3172
## r0[5] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 2719
## r0[6] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 2865
## r0[7] 0.00 0.00 0.01 0.01 -0.02 0.02 1.03 2699
## r0[8] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 2589
## r0[9] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 2598
## r0[10] 0.00 0.00 0.01 0.01 -0.02 0.02 1.03 2887
## r0[11] 0.00 0.00 0.01 0.01 -0.02 0.02 1.03 3730
## r0[12] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 2705
## r0[13] 0.00 0.00 0.01 0.01 -0.02 0.02 1.04 2441
## r0[14] 0.00 0.00 0.01 0.01 -0.02 0.02 1.04 3104
## r0[15] 0.00 0.00 0.01 0.01 -0.02 0.02 1.03 3611
## r0[16] 0.00 0.00 0.01 0.01 -0.02 0.02 1.03 3203
## r0[17] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 3015
## r0[18] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 3160
## r0[19] 0.00 0.00 0.01 0.01 -0.02 0.02 1.03 2525
## r0[20] 0.00 0.00 0.01 0.01 -0.02 0.02 1.04 2418
## m0[1] 3.19 3.19 0.02 0.01 3.17 3.21 1.00 1524
## m0[2] 3.26 3.25 0.02 0.01 3.23 3.28 1.00 2167
## m0[3] 3.56 3.56 0.02 0.02 3.53 3.59 1.00 1573
## m0[4] 3.32 3.32 0.01 0.01 3.30 3.34 1.00 2428
## m0[5] 3.31 3.31 0.01 0.01 3.28 3.33 1.00 2210
## m0[6] 3.43 3.43 0.02 0.02 3.40 3.45 1.00 2419
## m0[7] 3.30 3.30 0.02 0.01 3.28 3.33 1.01 2458
## m0[8] 3.11 3.11 0.02 0.02 3.09 3.14 1.01 1186
## m0[9] 3.40 3.40 0.02 0.01 3.37 3.42 1.00 2258
## m0[10] 2.82 2.82 0.03 0.02 2.77 2.86 1.01 858
## m0[11] 2.81 2.81 0.03 0.03 2.76 2.85 1.01 929
## m0[12] 3.46 3.46 0.02 0.02 3.43 3.49 1.00 1934
## m0[13] 3.46 3.46 0.02 0.02 3.43 3.49 1.00 1939
## m0[14] 2.98 2.98 0.02 0.02 2.95 3.01 1.01 995
## m0[15] 2.86 2.85 0.02 0.02 2.81 2.89 1.01 928
## m0[16] 3.28 3.28 0.02 0.01 3.25 3.30 1.00 2201
## m0[17] 3.15 3.15 0.02 0.02 3.12 3.17 1.01 1351
## m0[18] 3.51 3.51 0.02 0.02 3.49 3.54 1.00 2035
## m0[19] 3.50 3.50 0.02 0.02 3.47 3.53 1.00 2359
## m0[20] 2.86 2.86 0.02 0.02 2.82 2.90 1.00 901
## y0[1] 24.13 24.00 4.94 4.45 16.00 32.00 1.00 1976
## y0[2] 25.90 26.00 5.04 4.45 18.00 34.00 1.00 2003
## y0[3] 35.09 35.00 5.83 5.93 26.00 45.00 1.00 1867
## y0[4] 27.59 28.00 5.12 5.19 20.00 36.00 1.00 2074
## y0[5] 27.43 27.00 5.50 5.93 19.00 37.00 1.00 2060
## y0[6] 30.59 30.00 5.46 5.93 22.00 40.00 1.00 2094
## y0[7] 26.99 27.00 5.15 4.45 19.00 36.00 1.00 1969
## y0[8] 22.53 22.00 4.71 4.45 15.00 31.00 1.00 1866
## y0[9] 29.72 29.00 5.56 5.93 21.00 39.00 1.00 2026
## y0[10] 16.84 17.00 4.16 4.45 10.00 24.00 1.00 2103
## y0[11] 16.65 16.00 4.08 4.45 10.00 24.00 1.00 1837
## y0[12] 31.83 32.00 5.68 5.93 23.00 42.00 1.00 1998
## y0[13] 32.10 32.00 5.76 5.93 23.00 42.00 1.00 2100
## y0[14] 19.59 20.00 4.45 4.45 12.00 27.00 1.00 1967
## y0[15] 17.32 17.00 4.27 4.45 11.00 25.00 1.00 2009
## y0[16] 26.48 26.00 5.05 4.45 18.00 35.00 1.00 1944
## y0[17] 23.39 23.00 4.94 4.45 15.00 32.00 1.00 2131
## y0[18] 33.44 33.00 5.83 5.93 24.00 43.00 1.00 1985
## y0[19] 32.86 33.00 5.73 5.93 24.00 43.00 1.00 1867
## y0[20] 17.49 17.00 4.24 4.45 11.00 25.00 1.00 2022
## ess_tail
## 88
## 284
## 1029
## 102
## 1053
## 728
## 1005
## 797
## 867
## 1070
## 962
## 989
## 552
## 693
## 944
## 915
## 618
## 808
## 919
## 876
## 656
## 1137
## 702
## 845
## 1259
## 1313
## 1539
## 1202
## 1523
## 1358
## 1303
## 1442
## 1053
## 280
## 303
## 1358
## 1222
## 605
## 289
## 1373
## 1540
## 1467
## 1408
## 291
## 1955
## 1995
## 1768
## 1979
## 2032
## 1882
## 1877
## 1855
## 1954
## 1926
## 1829
## 2019
## 2055
## 2067
## 1908
## 1904
## 1801
## 1961
## 1955
## 1954
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=exp(m)),xy,col='black')
categorical dist. with parameters following Dirichlet dist.
h~Dir(a), h[0,1], sum(h)=1, a[0,1], sum(a)=1
y~Cat(h), y[0,1], sum(y)=1ex15-1.stan
data {
int N;
int K;
array[N] int<lower=1,upper=K> y;
}
parameters {
simplex[K] a;
simplex[K] h;
}
model {
h~dirichlet(a);
y~categorical(h);
}library(gtools)
n=20
h=rdirichlet(n,c(0.5,0.3,0.2))
summary(h)## V1 V2 V3
## Min. :0.008 Min. :0.000 Min. :0.000
## 1st Qu.:0.130 1st Qu.:0.010 1st Qu.:0.004
## Median :0.413 Median :0.144 Median :0.020
## Mean :0.449 Mean :0.345 Mean :0.206
## 3rd Qu.:0.777 3rd Qu.:0.751 3rd Qu.:0.318
## Max. :0.996 Max. :0.985 Max. :0.992c0=c(1,2,3)
for(i in 1:n) y[i]=sample(c0,1,prob=h[i,])
table(y) |> prop.table()## y
## 1 2 3
## 0.50 0.35 0.15data=list(N=n,K=ncol(h),y=y)
mdl=cmdstan_model('./ex15-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -24.9423
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 9 -20.4197 7.57839e-05 0.000195176 0.6424 0.6424 12
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -20.42
## a[1] 0.39
## a[2] 0.35
## a[3] 0.26
## h[1] 0.52
## h[2] 0.35
## h[3] 0.13mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -29.50 -29.16 1.45 1.23 -32.29 -27.75 1.00 953 1423
## a[1] 0.36 0.34 0.18 0.20 0.10 0.68 1.00 1777 1354
## a[2] 0.34 0.32 0.18 0.19 0.08 0.65 1.01 1369 1235
## a[3] 0.30 0.28 0.16 0.17 0.07 0.60 1.00 1606 1421
## h[1] 0.49 0.49 0.11 0.11 0.31 0.66 1.00 3231 1478
## h[2] 0.35 0.34 0.10 0.10 0.19 0.52 1.00 2680 1702
## h[3] 0.16 0.15 0.08 0.08 0.05 0.31 1.00 2122 1328
beta binomial dist.
y~betaB(n,a,b): y~B(n,p), p~beta(a,b)ex15-2.stan
data {
int N;
int n;
array[N] int y;
}
parameters {
real<lower=0> a;
real<lower=0> b;
}
model {
a~cauchy(0,5);
b~cauchy(0,5);
y~beta_binomial(n,a,b);
}
generated quantities{
int y1;
y1=beta_binomial_rng(n,a,b);
}n=30
a=3
b=4
p=rbeta(n,a,b)
n0=10
y=rbinom(n,n0,p)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 2.25 4.00 3.87 5.00 8.00par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,n=n0,y=y)
mdl=cmdstan_model('./ex15-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -204.579
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 -199.196 0.00183778 0.00380116 0.9979 0.9979 9
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -199.20
## a 3.25
## b 5.10
## y1 1.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -197.07 -196.79 1.02 0.76 -198.91 -196.10 1.00 449 782
## a 5.60 4.86 3.28 2.21 2.37 11.03 1.01 493 543
## b 8.88 7.71 5.19 3.50 3.60 17.53 1.01 504 550
## y1 3.89 4.00 2.00 1.48 1.00 7.00 1.00 1755 1684
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 